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1.0  INTRODUCTION 


This  is  the  Third  Annual  Technical  Report  under  Air  Force  Office  of 
Scientific  Research  Contract  F49620-80-C-0088,  "Fundamental  Properties  of 
Soils  for  Complex  Dynamic  Loadings".  The  report  covers  the  contractual 
period  1  August  1982  through  31  July  1983,  but  the  work  discussed  was 
mostly  accomplished  after  submission  of  the  Second  Annual  Technical  Report 
in  April,  1983  [Dass,  Merkle,  and  Bratton  (1983)]. 

The  FY  1983  modification  to  the  basic  contract  statement  of  work 
contains  three  tasks: 

E.l.e.  Response  of  a  Clay  and  Silt  to  Laboratory  Boundary  Conditions 

E.l.f.  Soil  Element  Model  (SEM)  Analysis  of  Laboratory  Test  Data 

E.l.g.  Theoretical  Development/Modification  of  Constitutive  Model 

The  first  two  tasks  depend  on  laboratory  test  data  not  yet  available,  and 

will  be  reported  separately.  The  third  task  reads  as  follows: 

Results  of  the  previous  work  will  establish  a  framework  within 
which  the  material  models  can  be  evaluated,  and  illustrate  the 
limitations  of  existing  models.  The  theoretical  development  of 
improved  constitutive  models  can  develop  along  one  of  two 
paths.  The  first  Involves  the  development  of  a  new  model.  The 
second  involves  modification  to  existing  models  to  include 
effects  not  currently  present.  The  selection  of  the  preferred 
technical  approach  will  be  made  and  preliminary  work  begun  on 
the  new  modeling  procedure.  The  work  begun  last  year  on  pore 
pressure,  rate  effects,  and  shear  behavior  will  be  continued 
and  new  aspects  of  soil  modeling  which  come  to  light  will  be 
reviewed.  Where  Inconsistencies  arise  the  emphasis  will  be 
placed  on  matching  insltu  behavior  as  opposed  to  laboratory 
behavior.  The  Initial  theoretical  development  work  will  be 
checked  In  an  ongoing  fashion  utilizing  the  Soil  Element  Model 
and  CIST  calculations. 

This  report  contains  three  major  sections.  The  first  section  deals 
with  the  general  equations  for  dynamic  response  of  a  saturated  soil,  which 
establish  the  mathematical  and  computational  framework  into  which  any  soil 


constitutive  model  must  fit.  The  second  section  deals  with  those  aspects 
of  soil  stress-strain  behavior  which  a  soil  constitutive  model  may  have  to 
reproduce.  The  third  and  final  section  presents  the  equations  of 
elastoplastlclty  needed  for  model  development,  and  then  discusses  the 
Initial  phase  of  that  development,  viz  the  shear  failure  criterion.  The 
proposed  shear  failure  criterion  has  several  convenient  features: 

1.  It  Is  related  to  stress  through  the  first  total  stress  Invariant 
and  the  second  and  third  deviator  stress  invariants,  each  of 
which  has  a  simple  physical  interpretation. 

2.  Its  parameters  can  be  determined  from  simple  linear  plots. 

3.  The  model  can  match  unequal  friction  angles  in  tri axial 
compression  and  extension. 

4.  The  ratio  of  octahedral  shear  to  octahedral  normal  stress  can  be 
calculated  directly  {without  iteration)  when  the  value  of  Lode's 
parameter  is  known. 

Other  features  of  the  model  are  yet  to  be  determined. 
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2.0  DYNAMIC  RESPONSE  OF  SATURATED  SOIL 
2.1  Effective  Stress 

In  a  saturated  soil  having  discrete  grains  with  negligible 
intergranular  contact  areas  each  grain  is  completely  surrounded  by  pore 
fluid.  Therefore  the  pore  pressure  acts  throughout  a  soil  element,  in 
both  the  pore  fluid  and  the  solid  grains.  Superimposed  on  the  hydrostatic 
intragranular  stress  acting  within  each  solid  grain  due  to  the  pore 
pressure  is  the  additional  intragranular  stress  due  to  intergranular 
forces.  A  plane  section  through  soil  will  generally  cut  through  both 
solid  grains  and  fluid-filled  pore  space.  The  intragranular  stress  acting 
on  the  solid  portion  of  the  section  must  balance  the  pore  pressure  acting 
on  the  cut  grain  surfaces  on  one  side  of  the  section,  plus  the 
intergranular  forces  acting  on  the  same  cut  grain  surfaces.  The  pore 
pressure  acting  on  one  side  of  the  fluid  portion  of  the  section  simply 
balances  the  pore  pressure  acting  on  the  fluid  portion  on  the  other  side 
of  the  section.  This  situation  is  shown  in  Figure  2.1.  Summing  forces  in 
the  vertical  direction  yields 
♦ 

IFy  +  =  (  IC1y  +  pAt)  -  (pAt  ♦  7AS)  =0  (2.1) 

so  that  the  intragranular  normal  stress  component,  7,  equals  the  sum  of 
normal  components  of  intergranular  forces  divided  by  the  area  of  solids. 


It  turns  out  to  be  more  convenient  to  normalize  the  sum  of  normal 
components  of  intergranular  forces  with  respect  to  At,  rather  than  with 
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respect  to  A  so  we  define  the  effective  normal  stress,  a,  by  the 

w 

equation 


(2.3) 


Comparison  of  Equations  ( 2.2)  and  (2.3)  shows  that 


=  t  —  a 
a  =  t-  a  =  v- 

*s  r 


(2.4) 


Returning  to  Equation  (2.1),  if  we  define  the  total  normal  stress,  o,  by 
the  equation 


oA^  =  pA^  +  oAs  =  (p  +  a)A^ 


(2.5) 


a  =  p  +  a 


so  that 


(2.6) 


O  =  O  -  P 


(2.7) 


Summing  forces  in  the  horizontal  direction  in  Figure  2.1  yields 


ZFh  =  IC1h  '  TAs  =  0 


(2.8) 


so  that  the  intragranular  shear  stress,  t  ,  Is 


(2.9) 


The  effective  shear  stress  Is  obtained  by  normalizing  the  sum  of 
tangential  components  of  Intergranular  forces  with  respect  to  At,  rather 
than  with  respect  to  A$,  so  that 
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and  therefore 


—  T 
T  *  T~T-n 


(2.11 


2.2  Grain  Compression  Due  to  Pore  Pressure 
If  the  bulk  modulus  of  solid  grains  Is  K$,  then  because  the  pore 
pressure,  p,  acts  throughout  each  grain  It  Is  apparent  that  each  grain 
undergoes  a  compressive  volumetric  strain  due  to  pore  pressure  of  amount 
aV 

*•-4. 


vs  p  *- 


(2.12 


2.3  Solid  Skeleton  Compression  Due  to  Pore  Pressure 
If  each  grain  of  the  soil  skeleton  undergoes  a  compressive  volumetric 
strain  due  to  pore  pressure  of  amount  p/K$,  then  the  entire  soil 
skeleton  will  undergo  the  same  compressive  volumetric  strain,  since  the 
skeleton  Is  composed  of  grains  In  contact.  (This  argument  assumes  no 
grain  slippage.)  Thus 

AVt 

_(_t)  =  p  (2 .13 

Vp  * 


Note  that  the  skeletal  compressive  volumetric  strain  defined  by 
Equation  2.13  Is  unrelated  to  effective  stress.  It  is  similar  In  nature 
to  displacements  of  a  framed  structure  due  to  temperature  change. 

2.4  Grain  Compression  Due  to  Effective  Stress 
There  are  two  components  of  Intragranular  stress:  p  and  .  Each 
causes  grain  compression.  The  component  of  grain  compressive  volumetric 
strain  due  to  Intergranular  forces  (effective  stress)  is 


The  Intragranular  shear  stresses,  (1  /  j),  cause  grain  distortion, 
and  therefore  skeletal  distortion.  However,  this  is  assumed  to  have  no 
effect  on  skeletal  volume,  and  is  therefore  viewed  as  part  of  the  overall 
skeletal  shear  response.  The  above  discussion  does  not  include  dilatancy, 
which  is  caused  by  relative  displacement  between  grains  due  to 
intergranular  slipping,  a  separate  mechanism. 

2.5  Darcy's  Law  and  Fluid  Drag 

If  the  volumetric  flow  rate  of  pore  fluid  through  a  soil  cross 
section  of  total  area  At  is  Q,  the  discharge  velocity,  w,  is  defined  by 
the  relation 


•  • 

Of  course,  the  actual  fluid  particle  velocity,  v,  is  larger  than  w,  since 
the  actual  fluid  particle  velocity  Is  the  ratio  of  volumetric  flow  rate  to 
flow  cross-sectional  area. 


-S-  =  - 

nAt  n 


(2.16 


Under  steady  flow  conditions  it  has  been  found  that  the  discharge  velocity 
is  related  to  the  pore  pressure  gradient  by  Darcy's  equation 


V  -k„p. 


ir\j 


(2.17 


where  k^j  Is  the  permeability  matrix.  Under  steady  flow  conditions 
inversion  of  Equations  (2.17)  yields 


P,,  -  -kf‘  ij 


(2.18) 


Now  consider  a  pore  fluid  caused  to  flow  under  a  given  pore  pressure 
gradient  through  a  soil  element  of  unit  volume,  shown  schematically  In 
Figure  2.2.  The  net  force  exerted  on  the  bounding  surface  of  the  soil 
skeleton  Is 

fsi  «  -(1  -  n)p,f  (2.19) 

and  Is  Independent  of  the  nature  of  the  flowing  fluid,  and  of  whether  the 
flow  is  steady  or  variable,  since  the  pore  pressure  gradient  is  assumed 
fixed.  The  remaining  portion  of  the  pressure  gradient  under  steady  flow 
Is  the  drag  force  exerted  on  the  body  of  the  soil  skeleton,  which  is 

fd1  *  ‘p,1  '  'si  *  ”p*1  4  «  -nP.f  *  nkjj  (2-^0) 

Under  unsteady  f low  conditions  the  fluid  drag  force  exerted  on  the  body  of 
the  soil  skeleton  is  assumed  to  still  be  given  by  the  same  expression 

fd1 .  pkjj  i}  tt.tn 

2.6  General  Equations 

Up  to  this  point  normal  stresses  have  been  considered  positive  In 
compression.  However,  for  reasons  of  notaMonal  and  computational 
convenience.  It  turns  out  to  be  easier  ler  normal  stresses 

positive  In  tension,  as  well  as  longitu.  ns  positive  In 

extension.  Since  pore  pressure  Is  compressive  ly  nature,  however,  pore 
pressure  will  continue  to  be  considered  positive  in  compression.  Thus  the 
relation  between  total  stress,  pore  pressure  and  effective  stress  Is 
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The  equation  of  motion  for  the  soil  skeleton  must  consider  effective 
stress  (considered  to  act  over  the  total  area  of  a  soil  element),  pore 
pressure  (acting  over  the  bounding  area  of  solids),  fluid  drag  force 
(computed  per  unit  total  volume  of  a  soil  element),  gravity  (acting  on  the 
solid  mass)  and  skeletal  acceleration.  Referring  to  the  soil  element 
shown  In  Figure  2.3,  the  skeleton  equation  of  motion  Is 

j  -  (1  -  n)p,j  ♦  nk^jwj  ♦  (1  -  n)p$gi  *  (1  -  n^iij  (2.23) 

where  g^  is  the  1th  component  of  gravity. 

The  equation  of  motion  for  the  fluid  phase  must  consider  the  pore 
pressure  (acting  over  the  bounding  flow  area),  the  skeletal  drag  force 
(computed  per  unit  total  volume  of  a  soil  element),  gravity  (acting  on  the 
fluid  mass),  and  the  fluid  acceleration.  Here  It  Is  Important  to 
recognize  that  the  discharge  velocity,  w.  Is  measured  relative  to  that  of 
the  soil  skeleton.  Again  referring  to  the  soil  element  shown  In 
Figure  2.3,  the  fluid  equation  of  motion  is 

-np.,  -  nkjJij  *  npf9,  =  nCf(S,  *  ^-1 
or 

-I’m  -  k(j  "j  *  “f3t  *  Vi  *  I  Vi 

Addition  of  Equations  (2.23)  and  (2.24),  setting 

(1  -  n)p$  ♦  npf  *  p 


(2.24) 


(2.25) 


(2.26) 


yields 


°tj,j  +  p9<  *  pu1  +  pfw1 


(2.27) 


The  soil  skeleton  stress-strain  equations,  in  incremental  form,  are 

^ij  =  °1jkldckl  (2,28) 

where  de^  is  the  matrix  of  incremental  skeletal  strains  associated 
with  effective  stress,  and  is  the  skeletal  elastoplastic 

incremental  stiffness  tensor. 

From  Equation  (2.13),  the  matrix  of  incremental  skeletal  strains  due 
to  pore  pressure  is 

^®1J  *  "  ^  6ij  (2*29) 

Assuming  there  are  no  causes  of  skeletal  strain  other  than  effective 
stress  and  pore  pressure,  it  follows  that  the  matrix  of  incremental 
skeletal  strains,  dc^j.  Is  given  by  the  expression 

d*ij  ■  d*?j  *  <S  (2-30> 

so  that 

dci j  *  dEij  "  **1J  (2,31) 

Substitution  of  Equations  (2.31)  and  (2.29)  into  Equation  (2.28)  yields 

*°1J  *  Dijkl(dckl  +  akl}  (2,32) 

The  pore  fluid  stress-strain  equation  is  simple  in  concept  because 
the  pore  fluid  Is  assumed  to  have  a  constant  bulk  modulus,  K^,  and  zero 
shear  modulus.  Not  quite  so  simple,  however,  is  evaluation  of  the  rate  of 
fluid  volumetric  strain  at  a  point.  Note  that  the  equations  in  this 
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section  concern  the  displacement  of  a  point  In  the  soil  skeleton  whose 
Initial  coordinates  are  specified,  and  the  pore  fluid  velocity  at  the  same 
displaced  skeletal  point.  The  analysis  Is  Lagranglan  with  respect  to  the 
soil  skeleton,  but  neither  Lagranglan  nor  Eulerian  with  respect  to  the 
pore  fluid.  This  Is  because  Darcy's  Law  applies  to  skeletal  and  pore 
fluid  elements  occupying  the  same  point  in  space  and  time.  Thus  we  track 
neither  a  particular  fluid  particle  nor  the  fluid  velocity  at  a  fixed 
point;  rather  we  track  the  fluid  velocity  at  a  moving  skeletal  point  whose 
Initial  position  Is  specified.  Pore  fluid  velocity  Is  thus  a  skeletal 
quantity,  like  skeletal  displacement,  because  It  Is  a  vector  tied  to  a 
point  In  the  skeleton  whose  Initial  position  Is  specified.  Returning  to 
the  problem  of  obtaining  the  pore  fluid  stress-strain  equation,  we  use  the 
equation  of  pore  fluid  mass  conservation  to  express  the  pore  fluid 
volumetric  strain  In  terms  of  quantities  already  defined.  The  pore  fluid 
mass  conservation  equation  Is 

(2.33 


Expansion  of  Equation  (2.33)  yields 


pf,1w1  +  pfw1,1 


s  -np.  -  np , 


il-J 

Mii 


(2.34) 


The  rate  of  Increase  of  porosity,  n,  is  the  sum  of  three  terms: 

(+)  the  rate  of  skeletal  volumetric  strain, 

(♦)  the  rate  of  solid  grain  volume  decrease  due  to  pore  pressure,  per 
unit  total  volume  (see  Equation  (2.12), (1  -  n)p/K$ 

(-)  the  rate  of  solid  grain  volume  increase  due  to  effective  stress, 
per  unit  total  volume  (see  Equation  (2.14), 


and  the  total  rate  of  change  of  pore  fluid  density  is 


1  Dpf  i 

Thus  Equation  (2.34)  can  be  written  in  the  form 


dwi,i  =  "dEii  *  (1  '  n)^  +  '  n 

Equations  (2.32)  and  (2.36)  can  be  written  in  the  form 

*®1J  ‘  “3^  dp  s  Dijkldekl 
and 


(2.35) 


(2.36) 


(2.37) 


ao11  #n  .  1  -  n.  .  j  .  / 

"  *7  dp  s  de" 


Finally,  the  skeletal  strain-displacement  equations  are 


(2.38) 
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.*• 7  \  .'•  ;.'■  .'.■-  .■■  .  *.”1  •_.  "i  ■ »  ■.*  ■>  ’jv  v  .■  .-_ 


! 


v\ 


o 

s' 

s1 


5n 


■tj  *  “j,i> 

where  Is  the  1th  component  of  skeletal  displacement. 

Summarizing  the  above  results,  we  have  [cf  Zienklewicz  and 
Bettess  (1982)]: 

Effective  Stress  Definition 


° 1j  B  °1j  +  p61j 


Soil  Skeleton  Equation  of  Motion 

•u.j  - (1  -  ">pm  *  i*;}  v  11  • 

Pore  Fluid  Equation  of  Motion 

"p,i  ■  kij  "j  *  °fst  *  »f“l  *  I  »f5i 

Total  Density 
(1  -  n)p5  ♦  nPf  =  p 

Soil  Element  Equation  of  Motion 


(1  -  n)ps‘u’j 


JU,J 


*  pgf  =  pUj  *  Pfwf 


Soil  Skeleton  Incremental  Stress-Strain  Equation 

<cij '  dp  ■  °ijkideki 

Pore  Fluid  Mass  Conservation  Equation 

’  (?7  *  TfI|dp  ■  deH  *  ""i.i 


(2.39) 


(2.22) 


(2.23) 


(2.25) 


(2.26) 


(2.27) 


(2.37) 


(2.38) 


v' 


•v 

4 


fcwl 

**k*'ir 

a 

v-l 


■1 


im 

•Vi 

•-VJ 

■•Vi 

m 


*.  .i 
& 


*• 
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Skeleton  Strain-Displacement  Equation 

e1j  *  7  (u1,j  +  “j.l*  (2’39 

Emphasis  In  this  research  effort  during  the  last  three  years  has  been 
on  Improving  current  models  for  the  skeletal  elastoplastic  Incremental 
stiffness  tensor,  However,  It  is  Important  to  keep  In  mind  how 

the  D  tensor  fits  Into  the  complete  set  of  equations  for  the  dynamic 
response  of  saturated  soils,  which  are  needed  to  solve  complex  dynamic 
problems  Involving  wave  propagation,  liquefaction,  and  spall. 

2.7  Undrained  Behavior 

As  a  particular  example  of  the  application  of  the  above  general 

equations  to  the  response  of  a  saturated  soil  having  compressible  solid 

grains,  consider  the  case  of  undrained  loading.  Assume  the  skeletal 

stress-strain  response  to  be  Incrementally  linear,  with  skeletal 

contralned  modulus  H^,  coefficient  of  lateral  stress  at  rest  "ok-  and 

bulk  modulus  K^.  For  undrained  loading  we  assume  no  flow  (w^  =  w^  = 

••  ## 

Wj  =  0),  and  also  neglect  Inertial  effects  (u^  «  0).  The  governing 
equations  then  reduce  to: 


°1j  *  °1j  *  p51j 


(2.22 


D1jkl(dekl  +  ^  ^1* 


(2.40 


dc 


11 


^11  ,n 

*7  A 


♦ 


V^ixp 

s 


(2.41 


The  skeletal  Incremental  stiffness  matrix,  Dfjid,  is: 


•»  N 


kl 

11 

22 

33 

12 

23 

31 

21 

32 

13 

Hk 

KokMk 

KokMk 

0 

0 

0 

0 

0 

0 

KokMk 

Mk 

KokMk 

0 

0 

0 

0 

0 

0 

KokMk 

KokMk 

Mk 

0 

0 

0 

0 

0 

0 

0 

0 

0 

2Gk 

0 

0 

0 

0 

0 

0 

0 

0 

0 

2Gk 

0 

0 

0 

0 

0 

0 

0 

0 

0 

O 

CM 

0 

0 

0 

0 

0 

0 

0 

0 

0 

2Gk 

0 

0 

0 

0 

0 

0 

0 

0 

0 

2Gk 

0 

0 

0 

0 

0 

0 

0 

0 

0 

2Gk 

where 


2Gk  .  Mk(l  -  Kok) 


1  *  2Kok  "ic 


Equation  (2.40)  yields 


D1ikl 


(dckl  +  lit  6kl} 


D11kldek)  A  D11kk  ^ 

- -j - ♦  “Bp —  dp 


■*r 


(2.42) 


(2.43) 


(2.44) 


7.  ™.  ^  ^  ■  .5  Kn  *• 


.  ..  ..  .  .  . 


*  Kkdcn  +  it  dP 


Substitution  of  Equation  (2.44)  into  Equation  (2.41)  yields 

dcn  *  ^  (Kkdeii  4 1^ dp)  •  4  nqJ1,dp 


or 


or 


where 


U  -  ^)de^  «  -[£-  (1  -  jq)  *  "(j^  -  ^)]dp 


de 


ii 


1  1 
1 

*i  \~*r 


1  - 


V  J 


dp.-  £ 


+  n 
*s 


t — r 

1 

tl 


1  JL 

lm*i 


so  that 


dp  *  -Fdc^j  *  -F«klde 


kl^kl 


Substitution  of  Equation  (2.22)  Into  Equation  (2.40)  yields 

D, 

■•s 


do1j  +  a1jdp  ■  Dijkldckl  4  dp 

Kt 


D1jkldtkl  4  1 r  41jdp 


or 


(2.45) 


(2.46) 


(2.47) 
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(2.48) 


Kk 

da1j  s  D1jkld£kl  '  (1  ‘  iq)4fjdp 

Substitution  of  Equation  (2.47)  into  Equation  (2.48)  yields 

Kk 

doij  =  Dijkldekl  +  (1  ‘  T^)Faijakldekl 

*  CDijkl  +  (1  "  l^)Fa1jak13dckl 
=  Dijkldekl  (2,49) 


where  the  undrained  elastoplastic  Incremental  stiffness  tensor, 
is  defined  as 


♦  (1 


Kk 

‘  lHFaijakl 


(2.50) 


For  the  hydrostatic  component  of  undrained  loading.  Equation  (2.49) 
yields 


+  (1 


■  »|,  *  <*  - 
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which  is  the  result  obtained  by  [Gassmann  (1951:  par.  59,  p.  15)]. 
For  constrained  (uniaxial)  compression.  Equation  (2.49)  yields 


«°i  -  V‘i  *  «  - 


=  +  (1  -  — ) F ]  de^ 


The  relation  between  effective  stress  and  total  strain  can  be 
obtained  by  substituting  Equation  (2.47)  into  Equation  (2.40),  or  directly 
from  Equation  (2.48): 


i  j  =  ijkl  kl 


VC.  30 


where 


*ijkl  =  Df jkl  ‘  iq  Fsij6kl 


(2.54 


Equations  (2.53)  and  (2.54)  can  be  used  to  construct  effective  stress 
paths  with  strain  contour  overlays. 

Equations  (2.51)  and  (2.52)  define  the  undrained  bulk  and  constrained 
moduli,  respectively. 


and  from  Equations  (2.55)  and  (2.56),  or  directly  from  Equation  (2.50)  we 
can  obtain  the  value  of  K0(J.  Equations  (2.55)  and  (2.56)  yield 

K.,  1  ♦  2K  . 


so  that 


'ou  11 


Equation  (2.50)  yields 


Mu-"k*  (‘  -  ^) F 
Vou  =  Vok  *  (*  •  it) 


so  that 


MkKok  +  (2  ‘  iq) 
Mk  +  (2  -  jr) F 


Combining  Equations  (2.47),  (2.51)  and  (2.55)  yields 


.  _$£_ 

-  AP,„  F  p 

a^77 ' 17^71  =  =  B 


The  Important  thing  to  notice  about  Equation  (2.57)  is  that  it  holds  even 
when  the  principal  total  stress  Increments  are  unequal.  Thus,  for  example 


A^2  *  /  AOj 


Equation  (2.57)  yields 


B  1 

=  -  yfiOj  +  2^0^ )  =  -  Ao^)] 


.v  V  V* 


3.0  SOIL  DYNAMIC  CONSTITUTIVE  MODEL  REQUIREMENTS 


3.1  The  Nature  of  Soil 

Soil  Is  a  particulate  material.  Soil  particles  vary  In  size,  shape, 
hardness  and  surface  texture,  and  although  they  can  be  bonded  together  by 
mineral  deposits,  this  Is  the  exception  rather  than  the  rule.  There  are 
four  primary  consequences  of  the  particulate  nature  of  soil  [Lambe  and 
Whitman  (1969:  Chapter  2)]: 

a)  Deformation  of  soil  Is  partly  the  result  of  Individual  particle 
deformation,  but  primarily  the  result  of  Interparticle  sliding 
and  rolling. 

b)  Soil  Is  Inherently  multiphase.  The  soil  particles  constitute  the 
solid  phase,  and  the  remaining  space  Is  pore  space.  The  pore 
space  Is  filled  by  pore  fluid,  consisting  of  a  gaseous  phase 
(usually  air)  and  a  liquid  phase  (usually  water).  In  dry  soil 
the  liquid  phase  Is  absent,  and  in  saturated  soil  the  gaseous 
phase  Is  absent.  The  pore  fluid  chemically  influences  the  nature 
of  soil  particle  surfaces.  Including  contact  surfaces,  and  hence 
affects  the  process  of  Interparticle  force  transmission  and 
resistance. 

c)  The  pore  fluid  can  flow  through  the  pore  space.  Whether  flowing 
or  still,  the  pore  fluid  physically  interacts  with  the  soil 
particles,  thus  further  Influencing  the  process  of  interparticle 
force  transmission  and  resistance. 

d)  Sudden  load  changes  are  carried  jointly  by  the  soil  skeleton  and 
the  pore  fluid.  The  resulting  change  In  pore  pressure  usually 


*n 


V-Vj. 
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causes  pore  fluid  flow,  which  alters  the  proportion  of  load 
carried  by  the  soil  skeleton  and  the  pore  fluid,  as  well  as 
changing  the  configuration  of  the  soil  skeleton. 

Because  soil  deforms  primarily  by  Interparticle  slip,  soil  strength 
Is  basically  frictional  in  nature;  and  because  pore  fluid  pressure  reduces 
Interparticle  contact  normal  forces,  the  strength  of  a  soil  element  Is 
controlled  by  the  difference  between  the  total  normal  stress  acting  on  the 
element  and  the  value  of  the  element  pore  pressure,  l.e.,  by  the  effective 
stress.  Because  of  the  nature  of  soil  formation  and  deposition  processes, 
natural  soil  deposits  are  often  inhomogeneous  and  inherently  anisotropic, 
and  soil  profiles  are  frequently  erratically  discontinuous. 

3.2  Soil  Stress-Strain  Characteristics 

Soil  stress-strain  characteristics  are  a  consequence  of  the 
particulate  nature  of  soil  and  the  processes  by  which  soils  are  formed, 
deposited  and  subsequently  altered  in  place.  The  following  list  of  soil 
stress-strain  characteristics  is  prioritized  for  construction  of  soil 
constitutive  models  to  predict  the  behavior  of  soil  masses  under  complex 
dynamic  loads,  such  as  explosions,  earthquakes,  and  moving  vehicles: 

a)  Soil  deformation  and  strength  are  governed  by  effective  stress. 

b)  Both  volumetric  and  deviatoric  stress-strain  curves  are 
nonlinear,  even  at  small  strains,  and  the  type  of  nonlinearity  is 
stress  path  dependent.  Figure  3.1  shows  the  continuous 
transition  from  concavity  to  convexity  with  respect  to  the 
vertical  strain  axis  of  a  plot  of  vertical  effective  stress 
versus  vertical  strain,  measured  in  a  drained  trlaxlal  test.  The 
parameter  controlling  the  shape  of  the  stress-strain  curve  is  the 
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direction  of  the  effective  stress  path.  At  mean  pressures  above 
500  psl  some  volumetric  stress-strain  curves  exhibit  a  convex 
yield  region  due  to  grain  crushing  at  highly  stressed 
Interparticle  contact  points,  but  at  even  higher  mean  pressure 
the  volumetric  stress-strain  curve  again  becomes  concave  to  the 
strain  axis.  Figure  3.2  Illustrates  the  above  behavior.  A 
similar  phenomenon  Is  observed  for  one-dimensional  compression 
curves  at  much  lower  stresses,  due  to  Interparticle  slip  followed 
by  subsequent  skeletal  stiffening. 

c)  Under  drained  conditions,  shear  strain  and  volumetric  strain  are 
coupled.  This  coupling  Is  called  dllatancy.  Under  undrained 
conditions  the  tendency  of  the  soil  skeleton  to  change  volume  is 
opposed  by  the  relative  incompressibility  of  the  pore  fluid, 
which  develops  an  excess  pore  pressure  sufficient  to  maintain  the 
soil  skeleton  at  constant,  or  near  constant  volume.  It  is  vital 
that  soil  dllatancy  be  correctly  modeled  in  order  to  obtain  the 
correct  pore  pressure  and  effective  stress  under  all  loading 
conditions. 

d)  At  large  shear  strain  a  given  soil  approaches  a  residual  or 
ultimate  shear  stress  and  void  ratio  which  depend  on  the 
confining  pressure,  but  are  Independent  of  the  Initial  void  ratio 
prior  to  shearing.  The  residual  or  ultimate  shear  stress  and 
void  ratio  define  the  critical  state  at  the  given  confining 
pressure  [Casagrande  (1936:  262)]. 

e)  In  approaching  the  critical  state  an  Initially  dense  or  over¬ 
consolidated  soil  will  attain  a  peak  shear  stress  greater  than 
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the  critical  state  value.  The  peak  stress  generally  corresponds 
closely  to  the  maximum  expansion  rate.  At  larger  shear  strains 
In  a  strain  controlled  test  the  shear  stress  decreases  (strain 
softening)  and  the  soil  continues  to  expand  at  a  decreasing  rate 
until  the  critical  state  Is  attained.  Both  dense  and  loose  soils 
show  a  tendency  toward  denslflcatlon  at  small  shear  strains,  due 
to  particle  rearrangement.  Loose  sands  Initially  compact,  then 
expand  as  they  approach  the  critical  state;  normally  consolidated 
clays  compact  throughout  their  approach  to  the  critical  state. 
Loose  sands  exhibit  steadily  increasing  shear  resistance  as  they 
approach  the  critical  state;  even  normally  consolidated  clays  can 
exhibit  a  peak  shear  resistance  with  subsequent  strain  softening 
as  they  approach  the  critical  state.  These  basic  soil  stress- 
strain  phenomena  are  Illustrated  In  Figures  3.3,  3.4,  3.5  and  3.6. 

f)  The  Intermediate  principal  effective  stress  can  have  a 
significant  Influence  on  both  the  peak  and  the  ultimate  friction 
angles.  Figures  3.7  through  3.26  [Merkle  (1971)]  show  soil 
strength  data  plotted  In  the  octahedral  plane.  In  those  plots  1> 
Is  the  Mohr-Coulomb  friction  angle,  and  u  Is  Lode's  parameter. 

If  og  had  no  Influence  on  ?,  the  data  points  would  all  lie  on  a 
straight  line  of  constant  ?.  More  will  be  said  about  these  plots 
In  Section  4. 

g)  Because  soil  particles  are  generally  not  bonded  together,  soil 
tensile  strength  Is  primarily  the  result  of  particle 
Interlocking,  and  is  very  small.  Soil  tensile  failure  causes 
stress  redistribution  In  a  loaded  soil  mass. 


h)  Plastic  (Irrecoverable)  volumetric  and  deviatorlc  strains  are 
both  generated  from  the  onset  of  loading. 

1)  Separate  yield  and  plastic  potential  functions  appear  to  be 
necessary  for  compression  and  shear,  for  a  classical  plasticity 
model.  Plastic  flow  Is  frequently  nonassoclatlve,  especially  In 
shear. 

j)  Soils  exhibit  the  Baushlnger  effect,  l.e.,  loading  beyond  the 
virgin  yield  point  In  one  direction  Increases  the  elastic  range 
and  yield  point  for  unloading  and  reloading  in  that  direction, 
but  decreases  the  elastic  range  and  yield  point  for  subsequent 
loading  In  the  opposite  direction  [Timoshenko  (1956  11:412); 

Nadal  (1950:20)]. 

k)  Soil  stress-strain  behavior  can  be  strain  rate  dependent,  both 
because  the  effective  stress-strain  behavior  of  the  soil  skeleton 
Is  strain  rate  dependent,  and  because  of  the  time  dependence  of 
pore  fluid  flow  and  the  associated  pore  pressure  adjustment. 

l)  Cyclic  loading  In  shear  and/or  compression  produces  a  number  of 
effects:  Initial  densiflcatlon;  hysteresis;  decreasing 
Increments  of  permanent  shear  and  volumetric  strain  with  each 
cycle,  leading  eventually  to  stable  hysteresis;  stiffening;  and 
decrease  In  damping. 

m)  Natural  soil  deposits  exhibit  both  Inherent  (depositional )  and 
stress-  (or  strain-)  induced  anisotropy. 

n)  Sample  disturbance  often  makes  the  stress-strain  behavior  of  a 
soil  sample  different  In  the  laboratory  from  what  it  would  have 
been  In-si tu. 
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The  Soil  Element  Model  has  been  and  will  continue  to  be  used  In  this 
research  program  to  test  the  ability  of  soil  constitutive  models  to 
reproduce  the  above  stress-strain  characteristics.  These  characteristics 
significantly  Influence  the  response  of  a  soil  mass  to  complex  dynamic 
loads  associated  with  explosions,  earthquakes  and  moving  vehicles. 


4.0  ELASTOPLASTIC  MODEL  DEVELOPMENT 

4.1  Basic  Equations 

If  a  cylindrical  soil  specimen  Is  consolidated  under  an  Isotropic 
stress  (<jjc  =  72c  =  ®3C)»  then  subjected  to  drained  compressive 
loading,  unloading,  and  reloading  under  constant  confining  stress 
(<r2  =  ^3  =  ®3C;  >_  ^3C)»  the  stress-strain  curve  appears  as 

shown  In  Figure  4.1.  Several  Important  features  are  shown  In  Figure  4.1: 

1.  The  stress-strain  curve  Is  nonlinear,  even  for  small  stresses  and 
strains. 

2.  Upon  unloading  from  point  A,  some  of  the  total  strain  is 
recoverable  (BE),  but  the  remainder  Is  irrecoverable  (OB). 

3.  Reloading  occurs  along  a  path  (BC)  somewhat  different  from  the 
unloading  path  (AB),  until  reaching  the  previous  maximum  stress. 
At  this  point  additional  loading  approaches  and  proceeds  along 
what  appears  to  be  a  continuation  of  the  virgin  compression  curve 
(0A),  with  little  apparent  Influence  of  previous  unloading  or 
reloading. 

4.  Unloading  and  reloading  occur  along  paths  whose  secant  from  zero 
to  maximum  stress  has  a  slope  very  close  to  that  of  the  initial 
tangent  to  the  stress-strain  curve.  This  means  that  the 
Irrecoverable  portion  of  any  strain  increment  is  essentially  the 
difference  between  the  total  strain  increment  and  the  strain 
Increment  associated  with  a  straight  line  stress-strain  curve 
having  a  slope  equal  to  that  of  the  initial  tangent  to  the  actual 
stress-strain  curve. 
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By  convention,  recoverable  strains  are  called  elastic,  and 
irrecoverable  strains  are  called  plastic.  If  the  unloading  and  reloading 
curves  in  Figure  4.1  both  retraced  the  virgin  loading  curve  (OA)  instead 
of  following  curves  (AB)  and  <BC )  the  stress-strain  behavior  would  be 
called  nonli nearly  elastic.  As  it  is,  the  linear  portion  of  the  stress- 
strain  behavior  shown  In  Figure  4.1  is  elastic,  and  the  nonlinear  portion 
is  plastic.  Since  some  of  the  stress-strain  behavior  is  elastic  and  the 
rest  Is  plastic,  the  overall  stress-strain  behavior  is  called 
elastoplastic. 

Multiaxial  elastoplasticity  theory  extrapolates  the  above  one- 
dimensional  stress-strain  observations,  and  assumes  that  plastic  strains 
are  superimposed  on  elastic  strains  calculated  according  to  the  theory  of 
elasticity.  Thus 


E1j  =  c1j  *  c1j 

where 


(4.1) 


(4.2) 


When  the  elastic  behavior  Is  isotropic.  Equation  (4.2)  reduces  to  the  form 

°ij  -  "Vkk'u  *  "i1  -  <4-3' 

where  M  ■  constrained  elastic  modulus 

Kq  «  coefficient  of  elastic  lateral  stress  at  rest. 

The  parameters  M  and  KQ  are  assumed  to  be  constant,  independent  of 
strain. 

From  Equation  (4.2)  it  follows  that 

do1J  *  D1jkldekl  (4,4) 
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so  that  a  unique  relation  exists  between  the  elastic  strain  Increments  at 
a  point  and  the  corresponding  stress  Increments.  However,  what  Is  known 
In  a  strain  controlled  formulation  Is  not  the  elastic  strain  Increments, 
but  the  total  strain  Increments,  dc.^,  which  differ  from  the 
elastic  strain  Increments  according  to  the  incremental  form  of 
Equation  (4.1), 

dE1j  *  de1j  +  dE1j  (4*5) 

Once  the  possibility  of  plastic  strains  is  recognized,  four  questions 
arise: 

1.  Can  plastic  strains  occur? 

2.  If  they  occur,  what  will  be  their  relative  algebraic  values? 

3.  If  they  occur,  what  will  be  their  actual  algebraic  values? 

4.  Will  they  occur? 

Obviously,  Question  2  is  a  subset  of  Question  3.  The  reasons  for  listing 
the  two  questions  separately  are  the  mathematical  and  physical  conditions 
used  to  answer  them,  which  are  explained  below. 

The  mathematical  theory  of  elastoplasticity  presented  here  contains 
four  parts,  each  needed  to  answer  one  of  the  above  four  questions: 

1.  A  yield  criterion,  assumed  to  be  of  the  form 

f{o1j*  c?j}  -  0  (4'6) 

satisfaction  of  which  Is  a  necessary  condition  for  the  occurrence 
of  additional  plastic  strain  at  a  point. 

2.  A  plastic  potential  function,  g(ojj),  for  which 

deP  «  dx  4s-  (4.7) 

1j  a°ij 


■m 


5 


V'.v 
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which  gives  the  relative  algebraic  values  of  the  plastic  strain 
Increments,  i.e.  the  direction  of  the  plastic  strain  increment 
vector  in  stress  space.  Equation  (4.7)  Is  called  a  flow  rule, 
and  the  scalar  constant  dx  is  determined  by  the  next  condition. 

3.  The  requirement  that  Equation  (4.6)  be  satisfied  not  only  at  the 
beginning  of  yielding,  but  throughout  yielding  as  well,  so  that 


df  =  don  ♦  dc?,  =  0 
3oij  acPj. 


(4.8) 


Equation  (4.8)  is  called  the  consistency  condition,  and  yields 
the  value  of  dx  in  Equation  (4.7).  It  therefore  permits 
calculation  of  the  actual  algebraic  values  of  the  plastic  strain 
increments. 

4.  The  requirement  that  the  calculated  plastic  strain  Increments 
lead  to  a  positive  plastic  work  increment, 

dWP  =  a^dejj  >  0  (4.9) 


If  Equation  (4.9)  Is  not  satisfied,  then  additional  plastic 
strain  does  not  occur  at  a  point,  in  which  case  all  strain 
increments  are  elastic. 

Equations  (4.6),  (4.7),  (4.8),  and  (4.9)  have  been  written  assuming 
one  yield  criterion  (or  yield  surface),  and  one  plastic  potential 
function.  There  can,  however,  be  more  than  one  yield  surface,  and  an 
equal  number  of  corresponding  plastic  potential  functions.  If  this 
happens,  then  the  above  four  equations  apply  to  each  active  yield 
surface.  Thus,  If  m  yield  surfaces  are  active,  there  will  be  a  set  of 
plastic  strain  Increments  for  each  active  yield  surface,  the  values  of 
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which  are  determined  by  4m  equations  (counting  Equation  (4.7)  as  one 
tensor  equation.) 


(4.10) 


(4.11) 


Each  stress  component,  can  be  expresed  as  a  function  of  the 
three  principal  stresses,  o^,  o2,  c-j,  and  the  nine  direction  cosines 
of  the  three  principal  stress  axes  with  respect  to  the  arbitrary  Cartesian 
axes  used  to  define  the  o^.  However,  If  a  unit  vector  pointing  in  the 
direction  of  the  1th  principal  stress  axis  is  e^,  then  because  the  three 
principal  stress  axes  are  orthogonal,  we  have 


ei  ‘  (4-12) 

Equation  (4.12)  represents  six  Independent  scalar  equations  involving  the 
nine  principal  stress  axis  direction  cosines.  Thus,  there  are  only  three 


Independent  principal  stress  axis  direction  cosines.  Let  them  be  <*j, 
<*2,  and  a^.  Then  we  can  write 


°1j=  °1j(or  °2*  °3*  °1’  °2*  a3] 


(4.13) 


If  a  material  is  Isotropic,  the  dependencies  of  the  yield  function. 


f.  In  Equation  (4.6)  and  of  the  plastic  potential  function,  g,  In 


Equation  (4.7)  on  the  principal  stresses  o^,  <r2,  and  °3  are 
Independent  of  the  orientation  of  the  principal  stress  axes  with  respect 


to  the  coordinate  axes.  This  Is  a  specific  application  of  the  principle 
of  material  frame  indifference  [Malvern  (1969:389)].  This  means  that  not 
only  do  a^,  a2 ,  and  a3  not  appear  in  either  equation,  but  also  the 
stress  functions  which  do  appear  are  insensitive  to  subscript 
interchanges,  l.e.,  they  are  symmetric  functions  of  c^,  a2  and  a3. 

The  total  stress  invariants  1^,  12,  and  I3  satisfy  the  required 
conditions  of  symmetry.  They  are: 


=  on  +  o22  +  o33 


/  all 

°12 

+ 

°22 

°23 

■f 

°33 

°31 

\a21 

°22 

°32 

°33 

a13 

all 

'll 

°12 

a13 

CVJ 

°22 

°23 

31 

a32 

a33 

Equation  (4.7)  gives  the  relative  values  of  the  plastic  strain 
increments,  from  which  the  relative  values  of  the  principal  plastic  strain 
increments  and  the  orientation  of  the  principal  plastic  strain  increment 
axes  can  be  determined.  Since  the  plastic  potential  function,  g,  is  a 
function  of  the  three  total  stress  invariants,  1^,  I2,  and  I3,  given 
by  Equations  (4.14),  (4.15),  and  (4.15),  we  have 

9  =  gdj,  I2.  I3)  (4.17 


so  that  Equation  (4.7)  can  be  written  in  the  form 


(4.18 


Now  Equations  (4.14)  and  (4.15)  yield 


1  0  *  s  j 


(4.19) 


"(a22  +  a33) 


-(o33  +  aU] 


(°11  +  °22^ 


s  aji  "  I1*1J  =  °1j  "  !1*1J 


(4.20) 


To  obtain  the  derivatives  al^/ao^.  we  notice  that  if  the  matrix 
of  cofactors  or  signed  minors  of  the  stress  matrix,  a  is  denoted 
Z,  then  since  according  to  Equation  (4.16)  the  determinant  of  o  is 
I3,  the  inverse  of  denoted  o"1,  is 

IT 

-1  - 


(4.21) 


Now  the  Laplace  expansion  for  the  determinant  1 3  is 


*3  =  1  °ij  Ei.j 


(4.22) 


and  direct  expansion  will  demonstrate  that 


dl3  -1  T 

3^  '  1 1  j  '  J3  °1  j 

T  T,-l  T  -1 
*  J3  °ij  s  I3°ij 


(4.23) 


Substitution  of  Equations  (4.19),  (4.20),  and  (4.23)  Into  Equation  (4.18) 
yields 
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de?j  s  dxC!f^  6ij  +  l?2{ofj  '  Vij}  +  Ifj  h  °ij] 


(4.24) 


Since  o_~l  has  the  same  principal  axes  as  does  £,  It  follows  from 
Equation  (4.24)  that  dep  also  has  the  same  principal  axes  as  does  £. 

This  condition  Is  a  consequence  of  the  assumption  of  material  isotropy, 
and  not  an  Independent  assumption. 

A  convenient  assumption  concerning  the  dependence  of  the  yield 
function,  f.  In  Equation  (4.6)  on  plastic  strain  Is  that  f  Is  a  function 
of  stress,  o^,  and  plastic  work,  Wp,  where  plastic  work  is  In  turn  a 
function  of  plastic  strain  [Malvern  (1969:367)]. 


f(o1j’  E1j)  =  ftoij*  wP(eij)] 


(4.25) 


Now  Equation  (4.9)  can  be  written  In  the  form 


jP  .  J^Ih.P  = 


pT  deij  B  °1jdeij 


(4.26) 


so  that 


3WP 

3c1j 


(4.27) 


The  application  of  the  above  equations  can  now  be  outlined. 

4.2  Stress  Control  Formulation 

When  stress  Increments  are  prescribed,  the  elastic  strain  increments, 
def*.  are  calculated  from  the  equation 


V  T  .  *  1  +  v 

r  Vu  e  aij 


(4.28) 
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where  E  *  Young's  elastic  modulus 
v  *  Poisson's  ratio 

Provided  Equation  (4.6)  is  satisfied,  so  that  yielding  can  occur,  we  write 
Equation  (4.8)  in  the  form 

df  .  do. .  *  dWp  .  0  { 4.29) 

Since  the  invariants  Ij,  Ig,  and  Ij  are  homogeneous  functions  of 
degree  1,  2,  and  3,  respectively,  Euler's  theorem  states  that 
Equation  (4.18)  yields 


dWp 


’ijdE1j  =  h  *  2  ffj  l2  +  3  Ifj  V  s  hdx  (4*30) 


where 


h  ■  h  * 2  %  '*  ‘ 3  Sj  '3 


(4.31) 


Equation  (4.30)  can  be  verified  by  direct  expansion  of  Equation  (4.24). 
Thus  Equation  (4.29)  can  be  written  in  the  form 


df  =  d a,.  *  dx ( h  =  0 

3oij  ij  3WP 


(4.32) 


and  therefore  the  flow  rule  proportionality  constant  is 
if 

1 

dX  : 


doij 


Tf 

"  8WP 

The  plastic  strain  Increments  are  calculated  from  Equation  (4.7), 

dep  ,  dx  -ML 

1J  30ij 


(4.33) 


(4.7) 


and  the  plastic  work  Increment  is  calculated  from  Equation  (4.30) 


dW*  »  hdx 


(4.34) 


4.3  Strain  Control  Formulation 

When  strain  Increments  are  prescribed,  the  elastic  strain  Increments 
are  not  Immediately  known.  Combining  Equations  (4.4),  (4.5),  and  (4.7) 
yields 
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Substitution  of  Equation  (4.41)  Into  Equation  (4.43)  yields 


'ijkl 


(D 


„  -iS_)(_*f-De 

Ijpq  3opq  a°rs  pskl 


) 


where 


(D 


,ep 


=  D 


Ijpq  3a 


_i£_)(4L  De 


'ijkl  "  ijkl 


£1 


3ar$  "rskl 


) 


_iL  De  SSL  -  h  1L 

3a  POrS  3a  ,,D 

pq  rs  3Wk 


(4.45) 


(4.46) 


4.4  Mixed  Boundary  Value  Formulation 

When  a  complementary  combination  of  stress  and  strain  Increments  is 
prescribed  (e.g.  as  In  constrained  compression).  Equation  (4.32)  can  be 
written  In  a  form  which  is  a  combination  of  Equations  (4.32)  and  (4.40). 


df  =  C-r^-  da..]H  ♦  [-£L  D*  ,(dEl.  -  dx  -r-S— )  ]  . 
3o ^ j  ij  do  3o j j  Ijkl  kl  30^1  de 


♦  dX  (h  IL.)  =  o 

3WP 


(4.47) 


where  the  symbol  and  [  ](je  mean  summation  only  over  indices 

of  prescribed  stress  or  strain  Increments.  The  flow  rule  proportionality 
constant  is  therefore 


dx 


‘  3f  , 

Nj  ij  J 

♦ 

do 

3f  r*e  j 

To^J  u1jklQekl 

de 

'  3f  De  _ 3j 

j  Ijkl  30| 

rj 

cl 

-  h 

de  3WP 

(4.48) 
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Equation  (4.48)  reduces  to  Equations  (4.33)  or  (4.41)  in  the  cases  of 
prescribed  stress  or  strain  increments,  respectively.  Provided  hdx  is 
positive,  the  plastic  strain  increments  are  calculated  from  Equation  (4.7) 

d‘?i  *  d»  iJJj  <4-71 

and  the  elastic  strain  Increments  corresponding  to  the  prescribed  total 
strain  Increments  are  calculated  from  Equation  (4.5) 

Cdeij  =  deij  ‘  deij]dE  (4,49) 

At  this  point  we  know  some  stress  Increments  but  not  the  corresponding 
elastic  strain  Increments,  and  the  complementary  elastic  strain  increments 
but  not  the  corresponding  stress  Increments.  To  calculate  the  unknown 
stress  and  elastic  strain  Increments  we  write  Equation  (4.42)  in  matrix 
form  as 


(«■*>  J 1 

& 

|  -12 

t  «>•>  i 

L  — 

►  * 

— 

1 - 

— 

CM 

rS 

jg 

ail 

i& 

t  tde®}  2 

(4.50) 


where  an  overbar  Indicates  a  matrix  of  prescribed  quantities,  and  the 
partitioned  matrices  and  are  square  and  symmetric. 

Expansion  of  Equation  (4.5)  yields 


{do}  j 

■4 

{3.'}  j 

{dee}  2 

(4.51) 

{do}  2 

-  fill 

(3«e)  j 

*  2?2 

{dee}  g 

(4.52) 

Equation  (4.52)  yields 


(4.53) 


{dee)2  -  D^*-1  |  tfc}  2  -  D|j  Qee}  jj 
and  substitution  of  Equation  (4.53)  Into  Equation  (4.51)  yields 
(d.)  j .  ojj  ra.e>  j <■  d'2o|j  -1  ra,>  2  -  o|j  <3.e)  jj 

4 

-  (Dji  -  dJ2522"1521)  {^e>  j  ♦  DjggJl*1  tta }  2  (4.54) 


The  remaining  unspecified  total  strain  increments  are  calculated  from 
Equation  (4.5), 

Cdeij  =  deij  +  deij3do  (4-5) 


The  elastoplastic  Incremental  compliance  tensor,  F®?^,  given 
by  Equation  (4.38),  and  the  elastoplastic  Incremental  stiffness  tensor, 
Dijkl*  9*ven  Equation  (4.46),  are  both  functions  of  stress  and 
plastic  work  only,  and  are  therefore  Independent  of  which  stress  and 
strain  Increments  are  prescribed.  Note  that  since 


D1jkldEkl 


D1jklFklpqdopq 


(4.55) 


It  must  be  that 


D1jklFklpq  =  41p°jq 


(4.56) 


so  that  DeP  and  FeP  are  the  Inverse  of  one  another. 


4.5  Computational  Format 

For  computing  purposes.  It  Is  convenient  to  write  many  of  the  above 
equations  in  matrix  form,  as  was  done  with  Equations  (4.50)  through 
(4.54) .  First  we  set 
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I 


I 

a 

\ 

,v 


and 


{o*} 


'  °11 ' 
a  22 
ff33 
a12 

<  °23  * 
a31 
a21 
°32 
‘  °13  , 


(4.57) 


{e*} 


'  r  ' 
E11 


e22 
c33 
c12 
e23 
C31 
C21 
c32 
L  Ei3  J 


(4.58) 


Equations  (4.57)  and  (4.58)  Indicate  that  the  stress  and  strain  spaces  are 
nine- dimensional.  However,  since 


°ij  *  °ji 

and 

eU  *  ej1 


(4.11) 

(4.59) 


all  stress  and  strain  points  are  restricted  to  six-dimensional  coincident 
subspaces.  We  therefore  set 

40 


iM 


.... . .  ■-  if  -  » i 


Equations  (4.60)  and  (4.61)  give  the  correct  expressions  for  work  and  for 
the  derivatives  of  the  stress  Invariants. 


The  reason  for  defining  the  six-dimensional  stress  and  strain  vectors  by 
Equations  (4.60)  and  (4.61)  Is  that  If  the  9x9  elastoplastlc  Incremental 
stiffness  matrix,  D®p,  In  the  equation 

(do*)  =  Dep  {de*}  (4.69) 

Is  symmetric,  then  the  stress  and  strain  vectors,  {a}  and  { e}  defined  by 
Equations  (4.60)  and  (4.61)  are  the  only  stress  and  strain  vectors  for 
which  the  6x6  elastoplastlc  incremental  stiffness  matrix,  £ep,  In  the 
equation 

{da}=Cep{d£}  (4.70) 


will  also  always  be  symmetric. 

In  matrix  form,  the  equations  for  a  complementary  combination  of 
prescribed  stress  and  strain  increments  are  as  follows.  Let 
(da)  *  column  matrix  of  unknown  stress  increments 
{ del  =  column  matrix  of  prescribed  stress  increments 
{ dy)  *  column  matrix  of  prescribed  strain  increments 
(da)  a  column  matrix  of  unknown  strain  increments 
Then  Equation  (4.47)  can  be  written  in  the  form 


df  =  }  {de}  ♦ 


[{da}  ^11 


I  { dy  }  -  dx  {|i} 

da 


♦  dx(h  il-)  =  0 


3W^ 


(4.71) 


so  that  the  flow  rule  proportionality  constant  is  therefore 


dx 


{T?} 
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{ del 
T" 


{ir} 

3a 


ce 

hi 


{  dy  } 


'  re 

-11  <7$  ' 


h  TfifP 


(4.72) 


■«  1 
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The  reduced  elastic  stiffness  matrix,  C^,  appearing  in  Equations 
(4.71)  and  (4.72)  contains  only  those  elements  associated  with  the 
prescribed  strain  increments,  {dy},  and  the  corresponding  (unknown)  stress 
increments,  {da}.  Provided  hdx  is  positive,  the  plastic  strain  increments 
are  calculated  from  the  matrix  form  of  Equation  (4.7), 

{depl  =  dx  (4?)  (4.73) 

30 

and  the  elastic  strain  increments  corresponding  to  the  prescribed  total 
strain  increments  are  calculated  from  Equation  (4.5) 

{dye}  =  {dy}  -  {dyP}  (4.74) 
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(4.81) 


cep  =  Ce 


ce  {|a}  (il}  ce 

—  3 o  3o  — 

T 


<£> 


4.6  Material  Behavior 

So  far  the  discussion  of  the  theory  of  elastoplasticity  has  been 
mathematical.  Material  behavior  comes  into  play  in  selecting  the  yield 
function,  f,  and  plastic  potential  function,  g,  on  the  basis  of  test 
data.  Frequently  the  mathematical  form  of  the  yield  function  is  based  on 
strength  or  failure  data,  and  a  hardening  function  is  added  to  convert  the 
failure  criterion  to  a  yield  function.  [Newmark  (1960:24)]  pointed  out 
that  the  definition  of  failure  should  be  as  precise  as  is  the  resulting 
failure  criterion.  From  a  mathematical  viewpoint,  once  the  yield  function 
and  plastic  potential  function  are  defined,  failure  occurs  when  the 
elastoplastic  incremental  stiffness  matrix,  Cep,  defined  by 
Equation  (4.81),  becomes  singular,  i.e.  when 

|Cep|=0  (4.82) 

The  theory  of  plasticity  was  first  developed  for  metals,  for  which 
the  yield  function  is  often  independent  of  the  hydrostatic  stress 
component.  (However,  ductility  of  metals  often  increases  with  increasing 
hydrostatic  stress.)  What  this  means  mathematically  is  that  in  three- 
dimensional  principal  stress  space  metallic  yield  functions  which  are 
Independent  of  the  hydrostatic  stress  component  are  right  cylinders,  with 
their  axis  along  the  line  oj  »  =  <73,  called  the  hydrostatic 

axis.  In  such  cases  interest  naturally  centers  on  the  shape  of  the 
intersection  of  the  yield  surface  with  a  plane  normal  to  the  hydrostatic 
axis,  called  a  deviator  or  octahedral  plane,  and  having  the  equation 


■„v. 
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(4.83) 


°1  +  °2  +  °3  s  !1  =  constant 

Figure  4.2  shows  the  hydrostatic  axis  and  an  octahedral  plane  in  principal 
stress  space. 

For  soils  the  hydrostatic  stress  component  definitely  influences  the 
failure  surface,  but  often  in  such  a  way  that  octahedral  cross-sections  of 
the  failure  surface  at  different  values  of  Ij  are  geometrically  similar, 
increasing  in  size  as  a  function  of  Ij.  When  this  happens  all  strength 
data  can  be  plotted  on  a  single  octahedral  plot  by  normalizing  the  data 
with  respect  to  the  1^  size  function.  This  is  how  Figures  3.7  through 
3.26  were  obtained.  For  example.  The  Drucker-Prager  failure  surface 
assumes  that  octahedral  cross-sections  are  circular,  with  the  radius  a 
linear  function  of  Ij  [Drucker  and  Prager  (1952:158)]. 

The  geometrical  justification  for  plotting  shear  strength  data  in  the 
octahedral  plane  using  coordinates 


(4.84) 


(4.85) 


as  shown  in  Figure  4.3  can  be  found  In  [Merkle  (1971:346)].  Using  the 
form  shown  in  Figure  4.4,  strength  data  from  many  different  investigations 
can  be  compared  on  a  common  basis. 

Although  the  Mohr-Coulomb  friction  angle,  0,  and  Lode's  parameter,  n, 
are  useful  for  plotting  strength  data  In  the  octahedral  plane,  the 
invariant  quantities 
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(4.86) 


Joct  =  7~ 


■J** 2  i  a  t; 

Toct  =  V"T  =  T  V(al  -  °2  + 


(o2  -  o3)2  ♦  (o3  -  Oj)‘ 


(4.87) 


cos  3u 


3 

T" 

7TJ77 

</> 


(4.88) 


J3  =  (ol  "  °oct)(o2  ‘  °oct)(o3  -  °oct) 


(4.89) 


are  more  useful  fn  constructing  curves  to  fit  the  data.  The  variables 

°oct’  Toct  an(J  °  are  s^own  *n  F*9ure  ^*3,  where  oj,  02,  and 
represent  principal  stresses. 

4.7  Drucker's  Equivalent  Stress  Function 

Although  his  principal  concern  was  developing  an  Invariant  shearing 
stress-strain  relation,  rather  than  defining  the  shape  of  the  octahedral 
cross-secion  of  a  failure  surface,  [Drucker  (1949:352)]  proposed  what 
amounts  to  a  formula  for  variation  of  Toct  with  u.  He  proposed  the 
following  expression  for  an  equivalent  shearing  stress, t  : 

'eq  ’  Vt(1  -  2-2w|/J2>1/6  <«•*» 


where  from  Equation  (4.88)  we  have 


J3  cos23y 
77  s  ~T7ir~ 
J2 


1  +  cos6u 

1775 — 


(4.91) 


Substitution  of  Equation  (4.91)  Into  Equation  (4.90)  yields 


_ 'eg 

Toct  =  T.  1  ♦  cos6u.l/6 
11 - B - > 


(4.92) 
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Equation  (4.90)  defines  a  smooth  curve  In  the  octahedral  plane  for  which 
lines  of  symmetry  occur  at  30  degree  Intervals  [cf  Hill  (1950:18)].  It 
causes  the  value  of  TQct  to  be  the  same  in  tri axial  extension  (u  =  60°) 
as  In  tri axial  compression  (u  *  120*),  and  as  Hill  explains,  such  an 
octahedral  cross-section  corresponds  to  an  Isotropic  material  which  does 
not  exhibit  a  Baushinger  effect.  Equation  (4.92)  can  be  written  in  terms 
of  Lode's  parameter,  4,  instead  of  the  octahedral  polar  angle,  u,  by  using 
the  relation  [Merkle  (1971:733)] 


2 

cos  3u 


2/  q  2 

6.75  4  =  Jt  (9  H  ’ 

(3  ♦  /)' 


(4.93) 


so  that  Equation  (4.92)  takes  the  form 


(4.94) 


When  u  =  *1,  Equation  (4.49)  yields 


(4.95) 


and  when  4=0,  Equation  (4.494)  yields 
Toct  =  Teq 


(4.96) 


4.8  Topping’s  Failure  Criterion 

[Topping  (1955:186)]  propose^  a  Mohr  circle  type  relation  in  the 
octahedral  plane,  to  account  for  the  possibility  that  the  octahedral  shear 
stress  values  at  failure  in  tri axial  compression  and  tri axial  extension 
may  be  different  at  the  same  octahedral  normal  stress.  His  equation  is 
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(4.97) 
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cos  3u 

Equation  (4.97)  can  also  be  written  in  the  form 

J. 


where 


or 


where 


Toct  *  *  *  B 


'3 

17? 


A  =  0.5(t:c  +  Te) 
B  =  1.3(tc  -  tc) 


Toct  -  A  -  C 


m(9  -l  ) 
(3  ♦  u*)3/* 


C  =  0.5(tc  -  re) 


(4.98) 

(4.99) 

(4.100) 

(4.101) 

(4.102) 


Note  that  Equation  (4.101)  permits  evaluation  of  the  constants  A  and  C  by 
a  straight  line  plot. 

4.9  Kirkpatrick’s  Failure  Criterion 

[Kirkpatrick  (1957)]  performed  both  conventional  triaxial  and 
thick-walled  cylinder  tests  on  Loch  Aline  sand,  for  the  purpose  of 
determining  the  shape  of  the  failure  surface  octahedral  cross-section. 

His  results  are  shown  in  Figure  3.10,  and  he  concluded  that  the 
Mohr-Coulomb  failure  criterion  was  a  good  fit  to  the  data.  In  fact, 
Kirkpatrick  felt  the  Mohr-Coulomb  criterion  fit  his  data  so  well  that  he 
decided  not  to  modify  the  axial  load  capability  of  his  thick  cylinder 
device  to  obtain  M  values  other  than  those  shown  in  Figure  3.10. 
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Both  Topping's  and  Kirkpatricks' s  results  were  cited  by  [Newmark 
(1960:28)]  as  examples  of  failure  criteria  having  roughly  triangular 
octahedral  cross- sections,  and  unequal  [octahedral  shear]  strengths  in 
trlaxlal  compression  and  extension. 

4.10  Coleman's  Failure  Criterion 

[Coleman  (1960:182)]  proposed  a  failure  criterion  based  on  an 
Invariant  formulation  of  the  Mohr-Coulomb  failure  criterion  for  a 
cohesionless  material.  From  Figure  4.3  we  have  (for  d  =  o). 


Js  sinu, 

*°oct  / 
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oct 
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so  that 


sin  d  = 
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When  u  =  +1 
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Therefore,  If  we  replace  sin  u2  by  -Jz/2t  and  cos  u2  by 
-l/2(cos  Sog)1^3,  the  resulting  expression  will  match  Equation  (4.103) 


for  m  =  *1.  The  resulting  expression  Is 


jr  ( - ^  4 

sin  0  = -  Ooct  - 

^ '  (vt)  ^,cos  3“21' 

j3Jz 

"  ZI1  ~  (J3\  1/3 

^  \r) 

Clearing  of  fractions  and  squaring  yields 


‘2I1  /J 3^ 

3J2  s  ~T  +  ( T ’) 


s1n20 


which  Is  Coleman's  expression.  The  mathematically  convenient  feature  of 
Equation  (4.105)  Is  that  the  Invariants  appear  separately,  and  In  a 
numerator. 

Equation  (4.103)  can  also  be  written  In  the  form 


Toct  v2(3  +  pi  )  sin  0 


and  when  u  *  1  (trl axial  extension)  Equation  (4.107)  yields 


4.11  Lomlze  and  Kryzhanov sky's  Failure  Criterion 
[Lomlze  and  Kryzhanovsky  (1967)]  performed  stress-controlled  true 
trlaxlal  tests  on  a  quartz  sand  from  the  Volga  region,  using  u  »  -1,  0  and 
1.  They  defined  strength  as  the  peak  value  of  TQCt  on  a  plot  of  T()Ct 
versus  roct,  at  constant  o0Ct  and  n,  and  used  the  following  Invariant 
expression  as  their  empirical  failure  criterion: 


where 


Du  = 


Now 


where 


n  -  J_  T°ct 
k  ~y/6  °oct 

a  =  1.73 

w*  x  260 


I3  =  o^ct  (6^  cos  3u  Dk3  -  9Dk2+  1) 


cos  3u 


m(9  -  M2) 

(3  *7')375 


(4.110) 

(4.111) 

(4.112) 

(4.113) 

(4.114) 


(4.115) 
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so  that 


3  6-^/T cos  3uD.  3  -  9DkZ  +  1 


27 

1 - ZTT 


Equation  (4.110)  can  therefore  be  written  In  the  form 
(6  \fT cos  3uDk3  -  9Dk2  ♦  l)1,73  -  f27jgg73  D, 


(4.116) 


or 
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(4.117) 


(6  -JT cos  3uDk3-  9Dk2+  l)1*73 

Dk - - m - 

Equations  (4.115)  and  (4.117)  can  be  used  to  find  the  value  of  Dk  for  a 
given  value  of  u,  by  iteration.  Figure  3.24  is  a  plot  of  the  results. 

For  a  more  thorough  review  of  investigations  of  the  effect  of  the 
intermediate  principal  stress  on  soil  shear  strength,  see  [Merkle 
(1971:Chapter  6)]. 

4.12  Modified  Lade  Model 

[Lade  (1972:137,  138)]  simplified  the  failure  criterion  of  Lomize  and 
Kryzhanovsky  by  deleting  the  factor  Dk  in  Equation  (4.110).  The  failure 
criterion  which  Lade  fit  to  true  triaxial  test  data  on  Monterey  No.  0  Sand 
was  written  in  the  form 

I3  -  k1I3  =  0  (4.118) 

Using  this  failure  criterion  as  a  basis.  Lade  developed  an  elastoplastic 
constitutive  model  having  one  associative  and  one  nonassociative  yield 
surface.  His  equations  are  given  below,  from  [Lade  and  Nelson  (1981)]. 

The  associative  compressive  yield  surface  and  plastic  potential  function 
are 

fc  *  f'c(o)  -  f"(Wc)  *  0  (4.119) 

where 

f^  «  Ij2  ♦  2I2  *  o2  +  o2  ♦  a2  (4.120) 

fc  *  «-121t 

and 
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9c  -  fc 


(4.122) 


p.  =  atmospheric  pressure 

Cl 

The  nonassocfative  expansive  yield  surface  and  plastic  potential  function 


where 


fp  *  y°>  -  y  v  ■  ° 


(4.123) 


fp  -  -  27»^>" 


(4.124) 


fp  =  Tij  at  failure 


(4.125) 


fp  -  ac_b“p(^)1/q.  (Q  >  0) 


(4.126) 


9p  =  Ij  -  [27  ♦ 


2 '17'  J13 


(4.127) 


The  octahedral  cross-section  of  Lade's  failure  surface  can  be  computed 


as  follows 


=  friction  angle  for  triaxial  compression 


1  ♦  sin  tr 


(N.  ♦  2)' 


B  .  1  - 


(4.128) 


(4.129) 


(4.130) 
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A 


cos  3u 


(4.131) 


’  f 

1 

2  s  Dk  *  =  I  ?B_-  Az  (iterate)  (4.132) 

oct 


Octahedral  cross-section  values  are  tabulated  for  6  =  30*  In  Table  4.1, 

v 

and  plotted  in  Figure  4.5.  This  plot  happens  to  be  an  excellent  fit  to 
Von  Karman's  data  in  Figure  3.7. 

Lade  and  Nelson  state  (p.504)  that  their  "yield  function  defines  the 
stress  levels  at  which  plastic  strain  increments  will  occur"  [emphasis 
added].  The  paper  is  silent  about  testing  a  given  total  strain  increment 
to  see  whether  it  will  cause  additional  plastic  strain,  and  in  this 
respect  the  model  appears  to  be  deficient.  However,  the  deficiency  is 
easy  to  correct,  and  the  correction  is  discussed  below. 

It  is  convenient  to  let  the  index  1  refer  to  quantities  related  to 
the  collapse  yield  surface,  and  the  index  2  refer  to  quantities  related  to 
the  expansive  yield  surface.  In  addition,  define  the  {s}  matrix  by 
Equation  (4.60),  the  {e }  matrix  by  Equation  (4.61),  and  the  6  x  6  elastic 
stiffness  matrix,  £e  based  on  Equation  (4.44),  i.e.. 
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(4.133) 
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where 


{de}  -  {dee}  ♦  {d£p} 


{deP}  =  {dePl }  ♦  {dep2} 


(4.134) 


(4.135) 


{do}  =  Ce  {dee} 


(4.136) 


Assuming  both  yield  surfaces  are  active  (which  may  not  be  the  case),  each 
of  the  two  consistency  conditions  takes  the  form 

af1  T  3f1  aWpj  T  D 1 

dfj  =  {aT*>  {d°}  +  {de  } 


3fi  T  3fj 

<d”!  *  "j  Jr 


=  0,  (j  =  1,  2;  no  sum)  (4.137) 


Now  let 


r1j  "  3^ 
39.- 


(4.138) 


(4.139) 


{ de p)  =  G  {dx} 


In  addition,  let 


Djj  =  h 


U  =  1  Z5J  #1j 


(no  sum) 


(4.140) 


(4.141) 


Equation  (4.137)  can  now  be  written  for  a  strain  controlled  condition  In 
the  form 


(4.142) 


{df}  -  FT  {do}  ♦  r°j  {dx} 

«  FTCe  {{de}  -  G  {dA}  }  +  rDj  {dx}  =  {0} 
and  Equation  (4.142)  can  be  solved  for  {dA}. 

{dA}  s  (FTCeG  -  rDj)_1  FTCe  {de}  (4.143) 

To  determine  whether  both  yield  surfaces  are  active,  we  examine  the 
plastic  work  Increments, 

dWp<J  =  hjdAj  (j«l,  2;  no  sum)  (4.144) 

Now  let 

hij  *  hia1j  (no  sum^  (4.145) 

Then  Equation  (4.144)  can  be  written  in  the  form 

{dwP}  =  rhj{dA}  =  {Thj  (FTCeG  -  rD_j  )_1  FTCe}{dE}  =  Q  {de}(4.146) 

where 

1  *  rh_,  (LVg  -  r0_,  )_1  iV  (4.147) 

The  Q  matrix  is  2  x  6,  which  means  that  Equation  (4.146)  requires  the 
total  strain  Increment  vector  {de}  to  have  positive  dot  products  with  both 
the  vectors  i_Qj  j  and  JJj  i  *n  order  for  both  yield  surfaces  to  be 
active.  This  requirement  Is  shown  graphically  In  Figure  4.6.  The  Q 
matrix  depends  only  on  the  current  stress  and  both  the  collapse  and 
expansive  plastic  work.  Thus  It  is  possible  to  tell  beforehand  whether  a 
given  total  strain  Increment  will  cause  both  yield  surfaces  to  be  active, 
when  the  current  stress  point  lies  on  the  Intersection  of  the  two  current 
yield  surfaces. 

-+  -►  -4- 

If  de  •  Qj  and  de  •  (?2  are  not  both  positive,  then  each  yield 
surface  must  be  examined  separately  to  see  whether  It  is  active  alone. 
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Thus  there  are  four  possibilities,  as  shown  below.  However,  it  is  not 


1 — _ _ 

Surface  2(P) 

— - — _ 

Active 

Inactive 

Surface  1  (C) 

Active 

CP 

C 

Inactive 

P 

E 

clear  from  a  comparison  of  Equations  (4.72)  and  (4.143)  that  the  outcome 

T  e 

jf  the  above  tests  will  be  unique,  unless  the  matrix  F  £  G  is  diagonal. 

Provided  both  yield  surfaces  are  active,  substitution  of 
Equations  (4.134),  (4.140)  and  (4.143)  into  Equation  (4.136)  yields 

(da)  =  Ce  {dee}  =  Ce  {{del  -  {deP> }  =  Ce  {  {d e}  -  G  {dx }  } 

=  Ce  H  -  G(FTCeG  -  rD_,  )_1  FTCe]  (de  }  =  Cep  {dE  }  (4.148) 

where  the  elastoplastic  Incremental  stiffness  matrix,  £®P,  is 

Cep  =  Ce  -  CeG(FTCeG  -  rD_,  )_1  FTCe  (4.149) 

Figures  4.7  through  4.11  show  stress-strain  curves  for  Antelope 
Valley  Sand  using  the  Lade  model,  computed  by  the  ARA  Soil  Element  Model 
(SEM)  program  using  model  parameters  given  by  [Lade  (1981)].  Figure  4.7 
shows  hydrostatic  compression  with  unloading  and  reloading.  Since  Lade's 
model  unloads  and  reloads  elastically,  it  does  not  show  hysteresis.  The 
elastic  moduli  are  Independent  of  strain,  so  the  two  unloading  lines  in 
Figure  4.7  are  parallel.  Figure  4.8  shows  uniaxial  compression  with 
unloading.  Figures  4.9,  4.10  and  4.11  show  curves  for  three  Isotropical ly 
consolidated  drained  trlaxlal  compression  tests,  at  constant  cell  pressure 
equal  to  the  consolidation  stress.  Figure  4.9  shows  plots  of  principal 
stress  difference  versus  axial  strain.  Since  the  three  samples  were  at 


the  same  void  ratio  prior  to  consolidation,  the  sample  subjected  to  the 
lowest  consolidation  stress  (14  psl)  behaved  as  a  dense  sand,  Increasing 
In  volume  for  axial  strains  larger  than  about  5.5  percent.  The  sample 
subjected  to  the  highest  consolidation  stress  (71  psi)  behaved  as  a  loose 
sand,  decreasing  in  volume  throughout  shear.  Figure  4.10  shows  volumetric 
strain  plotted  against  axial  strain,  and  Figure  4.11  shows  the  hydrostatic 
component  of  stress  plotted  against  volumetric  strain.  Figures  4.7  and 
4.11  would  be  identical  if  the  soil  were  linearly  elastic.  Obviously  it 
Is  not.  Comparison  of  these  two  figures  emphasizes  the  stress  path 
dependence  of  the  stress-strain  behavior  of  Antelope  Valley  Sand. 

4.13  Model  Development 

The  modified  Lade  model  discussed  above  is  a  partly  nonassociative, 
isotropic  hardening  elastoplastic  model,  with  both  yield  functions  and 
both  plastic  potential  functions  related  to  stress  through  the  total 
stress  invariants,  1^,  I2,  and  I-j.  Using  the  total  stress 
invariants  has  the  mathematical  advantage  that  differentiation  with 
respect  to  total  stress  is  straightforward.  However,  Lade’s  model  has  the 
disadvantages  that  It  cannot,  in  general,  achieve  an  exact  fit  to 
different  friction  angles  in  trlaxial  compression  and  extension,  and  the 
total  stress  invariants,  I2  and  Ij,  lack  a  simple  physical 
interpretation. 

In  contrast  to  the  above  situation,  octahedral  strength  plots  of  the 
type  shown  in  Figures  3.7  through  3.26  and  Figures  4.3  through  4.5  do  have 
a  simple  physical  Interpetatlon.  They  relate  t  and  y;  or  t  t  and  u;  or 
J2  and  J3,  where  J2  and  J3  are  the  second  and  third  deviator 
stress  Invariants,  arising  in  the  solution  of  the  principal  stress 
characteristic  equation.  Thus  the  question  naturally  arises  whether  it 
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might  be  more  convenient,  especially  from  the  standpoint  of  physical 
Interpretation,  to  make  the  yield  and  potential  functions  explicitly 
depend  on  Ij,  J2,  and  J3,  rather  than  on  Ij,  I2,  and  I3. 

First  of  all.  It  is  well  known  that  the  plastic  work  Increment, 
dWp,  can  be  expressed  as  the  sum  of  a  volumetric  term  and  a  deviatoric 


or  dlstortlonal  term.  This  is  done  by  defining  the  stress  deviator 
components,  s^,  and  plastic  strain  deviator  components,  e?j,  by 
the  equations 


EP 

1j  6ij 


(4.150) 

(4.151) 


The  expression  for  the  plastic  work  increment  can  now  be  written  in 
the  form 


°kk 
=  T 


sudefj 


(4.152) 


The  first  term  in  Equation  (4.152)  Is  the  volumetric  plastic  work 
increment;  the  second  term  Is  the  deviatoric  or  distortional  plastic  work 
Increment.  So  far,  so  good;  Equation  (4.152)  suggests  that  relating  f  and 
g  to  volumetric  and  dlstortlonal  Invariants  has  a  physical  basis. 

The  main  question,  then,  Is^whether  the  flow  rule  will  have  a 
convenient  mathematical  form  If  we  write 

g  =  gflj,  J2,  J3)  (4.153) 

In  place  of  Equation  (4.18),  we  can  write 
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and  by  analogy  with  Equations  (4.20)  and  (4.23). 

3J  n 


as 


kl 


=  slk  =  skl 


aJ3 

“h"  Skl 


and  from  Equation  (4.150) 

3Skl  1 

a^7  =  6ik4jl  “  7  6 k  1 6 i j 

Equation  (4.54)  can  now  be  written  in  the  form 

dEij  =  dx[lf^  6ij  +  (Uj){slk)Uik6jl  '  7  6kl5ij} 

+  (U^)(Skl,({1k4jl  "  7  4kl4ij)] 

“  dxC  41j  *  Uj  s1j  +  Uj  (Sij  •  T-  6ij)] 

Equation  (4.158)  yields  the  volumetric  plastic  strain  increment, 

dc?i  * dx(3  lf-}  =  dx(U  1 

11  3I1  3ooct 

so  that  the  deviatoric  plastic  strain  increments  are 


6ij 
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(4.160) 


=  dx[  I ^  S1j  +  1^S  ij  ‘  T~  a1j)] 

Equation  (4.160)  seems  workable  enough,  so  we  proceed  to  Investigate  the 
form  of  the  expression  for  the  plastic  work  Increment,  dWp. 

Equation  (4.152)  yields 

dwP  =  dx  C|f-  Ij  *  (2  Uj  J2  *  3  13- J3)]  (4.161) 

The  first  term  in  brackets  on  the  RHS  of  Equation  (4.161)  is  the 
volumetric  term;  the  last  two  terms  comprise  the  distortional  term. 

Equation  (4.31)  thus  has  two  alternate  forms: 


which  means  that 


(4.163) 


Thus,  there  appears  to  be  no  reason  for  not  using  J2  and  J3  instead  of 
I2  and  I3  in  the  formulation  of  the  yield  and  potential  functions,  and 
good  physical  justification  for  doing  so. 

The  proposed  expansive  yield  criterion  is  of  the  form  of 
Equation  (4.123),  i.e. 

fp  -  V'ocf  V  V  -  vv  *  0  ,4-IM) 

where 
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fn  -  Hr^Hl  -  e  cos  3u)( — —  ♦  b) 
P  P, 


V* 

p.  =  atmospheric  pressure 

a 


at  failure 


At  failure,  substitution  of  Equations  (4.165)  and  (4.166)  Into 


Equation  (4.164)  yields 


(^p^Hl  -  e  cos  3u)  =  — 
pa  pa 


*b  1  * b  ftf) 


(4.165) 


(4.166) 


(4.167) 


The  parameters  a,  b,  and  e  can  be  determined  from  a  series  of  triaxlal 
compression  and  extension  tests.  For  triaxlal  compression  («  =  120°), 
Equation  (4.167)  reduces  to 


(jjj  Br) 


where 


k  -  1  ~  e 
K1c  ~  a 

4c  ■  b'i-h1) 

For  triaxlal  extension  (u  =  60°),  Equation  (4.167)  reduces  to 


(4.168) 


(4.169) 


(4.170) 
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(rjjfg) 

■•■ft?) 


where 


°oct  1  +  e  r,  A  t  /°oct\  i  t,  x  l  /°oct\ 


1  ♦  e 


le  -  a 


(4.171) 


(4.172) 


(4.173) 


Plots  of  Equations  (4.168)  and  (4.171)  are  shown  in  Figure  4.12.  Such 
plots  are  often  referred  to  as  Southwell  plots  [Timoshenko  and  Gere 
(1961:191)].  Having  determined  the  parameters  klc,  k2c,  kle,  and 
k2e»  the  parameter  b  can  be  calculated  from  the  expressions 


k2c  k2e 

It —  =  t — 
*lc  Kle 


which  provides  a  consistency  check.  The  parameters  a  and  e  can  be 
calculated  from  Equations  (4.169)  and  (4.172),  written  in  the  form 

klca  ♦  e  -  1 
kle.  -  e  .  1 

so  that 


le  *lc 


•  -  1  -  kica  =  1 


2klc  kle  '  klc 
clc  +  kle  *  kle  +  klc 


(4.174) 


(4.175) 


(4.176) 
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Having  calculated  the  parameters  a  and  b  (and  e)  from  trlaxial 
compression  and  extension  tests,  the  accuracy  of  the  assumed  octahedral 
cross- sect  *  one  can  be  Investigated  by  a  series  of  true  trlaxial 
tests.  If,  In  Equation  (4.167)  we  set 


(4.177) 


then  Equation  (4.167)  can  be  written  In  the  form 


-  e  cos3u 


(4.178) 


which  Is  the  equation  for  an  ellipse  In  polar  coordinates. 

Equation  (4.178)  can  be  written  In  a  linear  form  to  obtain  the  octahedral 
eccentricity,  e,  as  a  consistency  check  on  the  previously  determined  value 
from  Equation  (4.176) 


«  1  -  e  cos3« 


(4.179) 


A  plot  of  Equation  (4.179)  Is  shown  In  Figure  4.13 

Octahedral  cross-section  data  for  the  case  (b  =  o;  6  =  32°; 
fle  *  35*)  are  tabulated  In  Table  4.2  and  plotted  In  Figure  4.13.  The 
calculation  sequence  used  to  obtain  the  values  shown  In  Table  4.2  and 


Figure  4.13  is  as  follows 

.  -Jl 

tan  u.  *  — 

<■  M 

zjz  sin  I 


(4.180) 


(4.181) 


Z-Jz  sin  0e 


♦  sin 


(4.182) 
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a  = 


(4.183) 


2z  z 
c  e 

2  +  za 
c  e 


e  = 


zc  -  2e 
zc  +  ze 


z  = 


1  -  e  cos3u. 


_  y/3  z  si 


nu. 


s*n  6  =  2  -  z  cos'*: 


(4.184) 


(4.185) 


(4.186) 


The  forms  of  the  compressive  yield  criterion  (cap),  and  the  plastic 
potential  functions  for  the  proposed  model  are  yet  to  be  determined,  and 
will  be  the  object  of  major  effort  during  FY84. 
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5.0  SUMMARY 

Recognizing  the  load-deformation  response  of  a  soil  mass  Is  governed 
by  effective  stress,  this  report  begins  by  developing  the  general 
equations  for  dynamic  response  of  a  saturated  soil  element.  Particular 
attention  is  paid  to  the  physical  significance  of  the  parameters  and  the 
equations.  These  equations  are  the  framework  into  which  a  soil 
constitutive  model  must  fit.  As  an  example  of  the  application  of  the 
general  equations,  the  problem  of  determining  incremental  total  stress 
moduli  for  static  undrained  (no  flow)  conditions  is  examined,  for  both 
Isotropic  (hydrostatic)  and  constrained  (one-dimensional)  compression. 

The  results  for  Isotropic  compression  agree  with  those  of  Gassmann.  The 
results  for  constrained  compression  should  be  useful  for  studies  of  soil 
liquefaction  under  explosive  or  earthquake  loading. 

Types  of  soil  stress- strain  behavior  observed  in  tests  are  examined, 
to  see  what  features  a  dynamic  soil  constitutive  model  should  possess. 
Included  are  effective  stress  dependence,  nonlinearity,  stress  path 
dependence,  dilatancy,  crlticial  state  behavior  at  large  shear  strains, 
peak  strength  behavior  (or  lack  of  It),  Influence  of  the  intermediate 
principal  stress  on  shear  stength,  low  tensile  strength,  inelasticity, 
nonassoclated  plastic  behavior,  the  Baushinger  effect,  rate  dependence, 
hysteresis,  decrease  of  damping  with  the  number  of  cycles  of  reversed 
loading,  anisotropy,  and  sample  disturbance. 

The  soil  constitutive  model  proposed  for  complex  dynamic  loading  is 
an  Isotropic,  strain  hardening  elastoplastlc  model.  The  basic  equations 
of  elastoplaticlty  are  developed  in  this  report  for  the  purpose  of 
emphasizing  their  logical  structure.  A  direct  (non-iterative)  solution  is 
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developed  for  the  case  of  mixed  boundary  conditions,  because  several  of 
the  most  conrnon  laboratory  soil  tests  are  of  this  type  (e.g.  strain- 
controlled  triaxial  compression  and  confined  or  uniaxial  compression),  and 
also  because  mixed  boundary  values  occur  frequently  in  (dynamic  finite 
difference  and  finite  element  computer  calculations  of  the  type  in  which 
the  proposed  model  will  be  used.  Criteria  for  distinguishing  between 
loading  and  unloading  are  carefully  examined,  because  of  the  oscillatory 
nature  of  soil  stress-strain  response  under  explosive  and  earthquake 
loading. 

The  proposed  model  shear  failure  criterion  has  the  following 
convenient  features: 

1.  It  Is  related  to  stress  through  the  first  total  stress  Invariant 
and  the  second  and  third  deviator  stress  Invariants,  each  of 
which  has  a  simple  physical  Interpretation. 

2.  Its  parameters  can  be  determined  from  simple  linear  plots. 

3.  The  model  can  match  unequal  friction  angles  In  triaxial 
compression  and  extension. 

4.  The  ratio  of  octahedral  shear  to  octahedral  normal  stress  can  be 
calculated  directly  (without  Iteration)  when  the  value  of  Lode's 
parameter  Is  known. 

Other  features  of  the  model  are  yet  to  be  determined. 

The  objectives  of  the  next  year's  effort  are  to  complete  the  proposed 
model  formulation,  and  to  demonstrate  Its  ability  to  reproduce  important 
aspects  of  observed  soil  stress-strain  behavior  discussed  in  Section  3. 
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TABLE  4.1 


LADE  FAILURE  SURFACE  OCTAHEDRAL  CROSS-SECTION 
FOR  tc  =  30  DEGREES 
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TABLE  4.2 


PROPOSED  FAILURE  SURFACE  OCTAHEDRAL  CROSS-SECTION 
FOR  (b  =  0;  <TC  =  32  DEGREES;  te  =  35  DEGREES) 


u2 

z 

sin 

-1.0 

120.000 

0.60679 

■ 

-0.9 

117.457 

0.60589 

-0.8 

114.791 

0.60304 

B 

-0.7 
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0 

90.000 

0.51938 

0.636 

0.1 

86.696 

0.50681 

0.633 

0.2 

83.413 

0.49526 

0.628 

0.3 

80.174 

0.48500 

0.622 

0.4 

76.996 

0.47620 

0.615 

0.5 

73.898 

0.46894 

0.608 

0.6 

70.893 

0.46321 

0.600 

0.7 

67.994 

0.45897 

0.593 

0.8 

65.209 

0.45610 

0.586 

0.9 
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Vertical  Effective  Stress  Mohr  Circle  Radius 


Axial  Strain 


Figure  3. 


3  Drained  Stress-Strain  Curves  for  Loose  and  Dense  Samp! 
of  the  Same  Sand,  Under  the  Same  Constant  Confining 
Pressure  [after  Lambe  and  Whitman  (1969:131)]. 


Deviator  Stress 


overconsol idated 


normally  consolidated 


Axial  Strain 


Axial  Strain 


Drained  Stress-Strain  Curves  for  Normally  Consolidated 
and  Overconsolidated  Samples  of  the  Same  Clay,  Under 
the  Same  Constant  Confining  Pressure  [after  Lambe  and 
Whitman  (1969:302,  312)1. 


Deviator  Stress 


30° 


Data  from  Kirkpatrick  (1957) 
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30° 


from  W°o d  (1958) 


Data  from  Sowa  (1963),  Skemption  and  Sowa 
(1963)  and  Wade  (1963) 


Figure  3.15.  Effect  of  a2  on  the  Strength  of  Remolded  Weald  Clay 


Figure  3.16.  Effect  of  o2  on  the  Strength  of  Remolded  Commercial 


Kaolinite,  o. 
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Figure  3.17.  Effect  of  o ^  on  the  Strength  of  Remolded  Commercial 
Kaollnlte,  oc  *  25  Psl. 
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Data  from  Shibata  and  Karube  (1965) 
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Data  from  Bell  (1965) 


Figure  3.21.  Effect  of  Og  on  the  Strength  of  Standard  Ottawa  Sand 
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a  V  i* 


Data  from  Green  and  Bishop  (1969) 
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L«d«  Failure  Surface  Octahedral  Cross-Section 
for  4C  •  30  Degrees. 
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Figure  4.6.  Condition  for  Both  Yield  Surfaces  to  be  Active. 
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Antelope  Valley  Sand 
(Isotropic  Compression) 


Figure  4.7.  Pressure  versus  Volumetric  Strain 


Axial  Strain 


Figure  4.8.  Total  Axial  Stress  versus  Axial  Strain 
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.13.  Linear  Plot  for  Determining  the 
Octahedral  Eccentricity. 
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Fl9ure  4.14  pH/s 
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